Spatial Domain Superresolution Reconstruction 

of Images 


by 

Kaushlendra Singh Sisodia 


"TH 

£9 }2oo])r^ 

Sul 83^ 

DEPARTMENT OF ELECTRICAL ENGINEERING 

INDIAN INSTITUTE OF TECHNOLOGY, KANPUR 

February, 2001 



f 



Spatial Domain Superresolution Reconstruction of 

Images 


A Thesis Submitted 

in Partial FuMUment of the Requirements 
for the Degree of 
Master of Technology 


by 

Kaushlendra Siagh Sisodia 



to the 

Department of Electrical Engineering 

Indian Institute of Technology, Kanpur 


February, 2001 





iiniiii 



Certificate 


This is to certify that the work contained in the thesis entitled “Spatial Domain 
Superresolution Reconstruction of Imaged’ , by Kaushlendra Singh Sisodia, has been 
carried out under our supervision and that this work has not been submitted elsewhere 
for the award of a degree. 


(Dr. K. S. Venkatesh) 

Deptt of Electrical Engineering, 
Indian Institute of Technology, 
Kanpur. 


(Dr. Sumana Gupta) 

Deptt of Electrical Engineering, 
Indian Institute of Technology, 
Kanpur. 



Acknowledgement 


I would like to thank to my thesis supervisors for their best guidance, advice, tolerance and 
continuous encouragement during the thesis work and more than anything else their faith 
in me. I will cherish my association with them for years to come. They have allowed me 
to learn and pursue various other interests of mine without any interference and I am veiy 
happy to have been their student. 

I am indebted to my parents whose inspiration and upbringing could bring me to the 
portals of my present achievement. I am thankful to my friends especially Rajiv Bajpai, 
Amitava Mukheijee, Ashwani and Vinay, whose encouragement and continuous interaction 
lead me to successfully completion of my thesis. 

Finally I am thankful to Almighty and my saviour Bhagvan, without his blessing and 
guidance nothing would have been possible. 



Abstract 


Superresolution (SR) reconstiuction for a linear space variant point spread function (LSV 
PSF) is highly computationally intensive and generally ill conditioned. In this thesis we 
propose two algorithms for spatial domain reconstruction of SR images. The first algo- 
rithm includes recursive SR reconstruction of an image from several undersampled de- 
graded frames, blurred by a space varying medium, with better reconstruction results than 
the existing method. The second algorithm includes a computationally efficient method for 
separable PSF taking the fact that 2-D shift matrix and 2-D downsampling matiix are al- 
ways separable and for getting a stable and unique solution of the ill-conditioned equations 
we use empirical regularized least square method. 



Synopsis 


The supeiresolution(SR) reconstruction of a high resolution image is an image restoration 
method from several degraded low resolution frames The supeixesolution concept evolved 
because large telescopes have reached a limit of practical realization, restricted by the high 
costs of manufacturing large and precise optics. Superresolution overcame these limitations 
by processing low resolution images and reconstructing a high resolution image. The con- 
cept behind SR reconstruction is to extract different(extra) information from each degraded 
frame of low resolution frame and to use it to reconstruct a high resolution image. 

Superresolution reconstruction for a linear space variant point spread function (LSV 
PSF) is highly computationally intensive and generally highly ill conditioned. The results 
in the thesis, motivated by the need to satisfactorily process different satellite pictures, 
are, however applicable in other situations which are faced with space variant blurring and 
noise addition. The algorithms and techniques outlined provide noise reduction and an 
increase in spatial resolution using several blurred images of a common object. In this 
thesis we propose two algorithms for spatial domain reconstruction of SR images. The first 
algorithm includes recursive SR reconstruction of an image from several undersampled 
degraded frames, blurred by a space varying medium, with better reconstruction results 
than the existing method. The algorithm may be used for reconstmction of an SR image 
from several frames obtained from multiple cameras with different sensor element sizes 
with each frame degraded by a different PSF. Computational complexity can be further 
reduced by using iterative methods. 

The second algorithm includes a computationally efficient method for separable PSF 
taking the fact that 2-D shift matrix and 2-D downsampling matrix are always separable and 
for getting a stable and unique solution of the ill-conditioned equations we use empirical 
regularized least square method. 

After each update of the estimate, we use the a priori information that each pixel value 
in the source image is nonnegative. 
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Chapter 1 
Introduction 


1.1 Introduction and motivation 

Superresolution reconstruction is a method of image restoration. The field of image restora- 
tion began primarily with the efforts of scientists involved in the space progr ams of both 
the United States and the former Soviet Union in the 1950s and early 1960s. These pro- 
grams were responsible for producing many incredible images of the earth and our solar 
system that were hitherto unimaginable. Such images held untold scientific benefits which 
only became clear in ensuing years as the race for the moon began to consume more and 
more of scientific effort and budgets. However, the images obtained from the various plan- 
etary missions of the time, such as the Ranger, Lunar Orbiter, and Mariner missions, were 
subject to much photographic degradation. This was the result of substandard imaging 
environments, the vibration in machinery and the spinning and tumbling of the space craft. 

The superresolution concept evolved because large telescopes have reached a limit of 
practical realization, restricted by the high costs of manufacturing large and precise op- 
tics. Superresolution overcame these limitations by processing low resolution images and 
reconstmcting a high resolution image. 

Several Earth Resources Technology (ERTS) Satellites (renamed later as LANDSAT) 
were launched to obtain multispectral images for subsequent processing in order to identify 
different environment phenomenon such as the distribution and general type of vegetation, 
regional geological stmctures, and the area extent of surface water. Each LANDSAT flew 
in a circular orbit, about 570 noiles above the surface of earth, executing circles approxi- 
mately 14 times per day. Each daytime orbital pass was from north to south with repetitive 
coverage every 18 days or so. 

The superresolution (SR) reconstruction principle as applied upon several frames came 
somewhat later, with the launch of the widely publicized Hubble Space Telescope. Launched 
in 1990 as a joint project of the North American and European space agencies NASA 
and ESA, this optical observatory was designed to yield a breakthrough in astronomy by 
providing space images with unprecedented spatial resolution. Unfortunately, soon after 
the launch, engineers discovered a manufacturing defect in the main mirror of the tele- 
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scope, leading to severe distortion that could eventually be fixed only in late 1993. The 
astronomers improved the blurred images obtained from the Hubble craft by numerical 
reconstruction, i.e., by solving an inverse problem. 

The results of the thesis, motivated by the need to satisfactorily process different satel- 
lite pictures, are, however applicable in other situations which are faced with space variant 
blurring and noise addition. The algorithms and techniques outlined provide noise reduc- 
tion and an increase in spatial resolution using several blurred images of a common object. 

With the availability of frame grabbers capable of acquiring multiple frames of video, 
there is a growing interest in superresolution image reconstruction from video, as well as 
video reconstruction, whereby multiple frames are used to overcome the inherent reso- 
lution limitations of a low-resolution camera system. SR reconstruction proves useful in 
many practical applications, including printing SR stills from video, where it is desired to 
enlarge an image and increase the detail. Because video signals are commonly interlaced, 
creating SR stills requires a combination of deinterlacing and the removal of acquisition 
degradation. 


1.2 Superresolution Reconstruction 

The SR concept in image restoration as originally formulated around the 1970’s has now 
undergone an evolution, and the present-day idea of SR reconstruction is somewhat differ- 
ent even though the basic concept remains the same: 

An incoherent object scene 0{x) is a spatial radiance distribution; it therefore cannot 
be negative; 

0{x) > 0, all X. (1.1) 

Let us consider the problem of restoring 0{x) from knowledge of its degraded image /(y) 
and the point-impulse response S{x) characteristic of the imagery. 

Let the unknown object be an incoherent, planar radiance distribution. The object plane is 
imagined to be subdivided into J equal-sized, elemental cells of extent Az, and Xj locates 
the center of the jth cell. (We use one-dimensional notation, for simplicity.) It is logical to 
make Ax just the required limit of resolution desired by the user in the final restoration. It 
therefore suffices to represent the object as a sequence of average radiance values 0{xj) = 
Oj, j=l,2, . . . , J, over the elemental cells. The light from an object is assumed to be 
imperfectly imaged by a linear system such as a lens. An image plane intercepts the light 
and the formed image is sampled by a sensor array at points y(rn); m = 1, 2, 3, . . . , M. If 
the image is bandlimited, it obeys a sampling theorem. In this case it is natural to space the 
sensor points y(m) uniformly at the Nyquist interval. 

In order to permit a higher resolution in the restoration of 0{j) than in the image, 
object-cell size Ax must be smaller than the image-data interval. 

Ax < y{m) — y{m — 1) (1.2) 
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This is the basic concept of superresolution. We can present the superresolution concept in 
either of the following two ways: 

1.2.1 The Classical Concept 

The original SR problem was posed as follows: can photographic images such as those 
of star clusters be superresolved? Can already good images be improved by a restoring 
method? To study the question, objects have been prepared that can be resolved only if 
the bandwidth in the restoration exceeds that of image data. Hence superresolution was 
known as extrapolation of the signal, that could be applied to a space limited function (i.e. 
f{x) = 0 for \x\ > a) whose Fourier transform is given over a finite frequency band. Ex- 
trapolation of the spectrum of an object beyond the known frequency band of the imaging 
system is called Superresolution. The detected image I(y) of an incoherent object 0(x) is 
formed through a linear relation: 

I{y) = f 0{x)S(y - x)dx -f N{y) (1.3) 

J scene 

where x, y are space coordinates and O, I have imits of radiance and irradiance, respectively. 
N(y) is random noise that is characteristic of the particular method chosen for detecting 
the image. S(y) is the diffraction image of a point object, and is called the point spread 
function(psf). 

Now the problem is inversion of (1.3) for restoration of object 0(x), i.e., estimating 
0(x) when I(y) and S(y) are known at a given subdivision of points y, and to a given level 
of accuracy. We call the estimate of 0(x) its restoration, and denote it as 0(x). 

Two mam obstacles to perfect restoration of 0(x) are the presence of noise N(y) in the 
image data, and a sharp cutoff in the spectmm s(a;) of the PSF. If we denote the cutoff 
frequency for s(a;) as , Fourier inversion of (1.3) yields: 

i(aj) = o(ui)s(uj) -I- n(u) for jw] < (1.4) 

n(uj) for |a;| > Q. (1.5) 


Lower case quantities in (1.4) and (1.5) correspond to the Fourier transform of the capital- 
ized quantities in (1.3). Evidently, the existence of a finite causes the image spectrum 
to lack components at all frequencies |a;| > This seems to apply that these higher - 
frequency components of the object cannot ever be recovered in the restored image. If so, 
the resolution in the restoration would be limited to a smallest separation of about 

2*pi/Q = R, (1.6) 

depending on the choice of resolution criterion that is used. R is called the Rayleigh reso- 
lution length. 
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1.2.2 The Modem View 


During 1990’s superresolution idea became somewhat changed now instead of single frame, 
many frames were used for reconstruction of a high resolution image. Now the problem is 
to restore a high resolution image from a sequence of low resolution, blurred, undersam- 
pled, discrete frames, where each frame has been shifted with respect to a reference frame. 
This problem is highly ill posed. Due to which the the solution of inverse problem may not 
exist, unique, or continuous or any combination of these three cases. Hence for making the 
solution possible we use regularized methods. Concept of regularized method is discussed 
widely in further chapter. 

In terms of a degree of freedom(DOF) analysis we have to reconstmct an image with a 
higher degree of freedom equal to from undersampled degraded frames having a lower 
degree of freedom. The concept of the degree of freedom in an imaging system is basic to 
all signal processing and evolves from the desire to minimize the number of independent 
parameters that determine an image. Basically, the intuitive meaning of degree of freedom 
might be related to the concept of the number of independent samples or data points nec- 
essary to define an image. Thus, simply oversampling an image by a factor of 10 does not 
imply that we obtain an order of magnitude increase in degrees of freedom or independent 
data points, as the oversampled image pixels would have high mutual correlation. Thus, in 
badly blurred systems, one would expect a lower number of DOF than in a well-corrected 
system even though the image sampling rate remains fixed. Ideally, the DOF describing an 
imaging system would be an invariant measure of the maximum resolution achievable by 
such devices. The DOF measxirement would provide a quantitative evaluation of the num- 
ber of independent parameters or variables obtainable from a given system without concern 
to Nyquist rate and oversampling phenomena. 

Hence in the multiframes superresolution problem, a higher DOF image is (re)constmcted 
from lower DOF frames. In the recursive SR algorithm, the DOF of the reconstructed high 
resolution image increases after adding each new frame. 


1.3 Related works 

The superresolution restoration idea was first presented by Tsay and Huang [11]. They 
used the frequency domain approach to demonstrate the ability to reconstruct one improved 
resolution image from several downsampled noisefree versions of it, based on the spatial 
aliasing effect In [10] single image restoration algorithm is proposed using iterative con- 
strained least square estimation, which uses properties of a block circulant matrix; which 
could not be preserved for the case of undersampled frames. 

Using the aliasing relationship between the undersampled frames and the reference im- 
age, a weighted recursive least square theory based algorithm is developed in the wavenum- 
ber domain [12]. They considered frames noisy and undersampled but not blurred. A 
frequency domain recursive algorithm for the restoration of superresolution images from 
noisy and blurred measurements is suggested in [13]. This algorithm could be used only 
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for LSI case, where the blur, the motion, and the decimation are all space variant. 

The work on multi-channel superresolution has been done in [14] which is a mean of 
recovering high frequency information by trading off the temporal bandwidth. Since the 
problem is usually ill-posed, one needs to impose some regularity constraints. However, 
regularity constraints tend to attenuate the high frequency contents of the data (usually 
present in the form of discontinuities). They investigated this inherent contradiction be- 
tween regularization and super-resolution issue in the context of adaptive regularization, 
using functions (convex, non-convex, bounded, unbounded). 

Conjugate gradient methods for superresolution has been shown to accelerate conver- 
gence to the solution [15]. They utilized a combination of Tikhonov-Miller regulariza- 
tion and positivity constraints as a means of regularizing the conjugate gradient algorithm. 
R.R.Schultz and R.L.Stevenson proposed a Bayesian motion estimation technique which 
modeled the motion field with a discontinuity-preserving prior [16]. 

In [17] a unified methodology is proposed utilizing the previous image restoration tools, 
the maximum likelihood (ML) estimator, the maximum a posteriori probability (MAP) 
estimator, and the set theoretic approach using projection onto convex sets (POCS). 

Michael Elad and Arie Feuer introduced the problem of reconstructing a super-resolution 
image sequence from a given low resolution sequence [18]. They proposed two iterative 
algorithms, the R-SD and the R-LMS, to generate the desired image sequence. These algo- 
rithms assumed the knowledge of the blur, the down-sampling, the sequences motion, and 
the measurements noise characteristics, and apply a sequential reconstruction process. 


1.4 Application of Superresolution Reconstruction 

Single image restoration has become a classical chapter in image processing theory, with a 
direct generalization to the restoration of continuous image sequences.The term continuous 
corresponds to the basic assumption that the image sequence contains one filmed scene. A 
standard video camera can be viewed as a source for such signals. Suppressing an additive 
noise in a continuous image sequence is an important preprocessing stage in many appli- 
cations such as image sequence coding and computer vision algorithms. Deblurring such 
a signal is an important tool for the enhancement of visual data presentation to the human 
viewer. 

In many applications such as in astronomy (satellite remote sensing) and computer vi- 
sion applications, it is required to reconstruct a high-resolution image from multiframes of 
undersampled low resolution, blurred and noisy images. Hence the scope and application 
of such superresolution reconstmction methods arises in following areas. 

(1) Remote sensing: where several images of the same area are given, and an improved 
resolution image is sought. 

(2) Frame freeze in video: whereas a typical single frame in the video signal is generally 
of poor quality and is not suitable for hard copy printout. Enhancement of a freeze image 
can be effected by using several successive images merged together by a superresolution 
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algorithm. 

(3) Medical imaging (CT, MRl, ultrasound, etc.}: these involve the acquisition of several 
images, that are limited in resolution quality. In the area of quantitative autoradiography 
(QAR), images are obtained by exposing X-ray sensitive film to a radioactive specimen. 
QAR is performed in post-mortem studies, and provides a higher resolution than tech- 
niques such as positron emissiontomography (PET), X-ray computed tomography (CAT), 
and magnetic resonance imaging (MRI). 

(4) Conversion from NTSC video to HDTV standard: printing from a NTSC source and 
conversion of NTSC source material to high definition television (HDTV) format are some 
recent application that motivate SR image and video reconstruction from low resolution 
(LR) and possibly blurred sources. 


1.5 Organization of the thesis 

hi Chapter 1 we discuss the motivation and application of the work, different conceptions 
of superresolution and review the work done previously in this area. In Chapter 2, we ex- 
plain the problem formulation and our algorithm. In Chapter 3, the inverse problem and 
its solution by least squares regularization are discussed. In Chapter 4 we discuss adap- 
tive Wiener filtering and recursive reconstruction procedure. In Chapter 5, we discuss the 
results. Finally, in last in Chapter 6, we conclude the work and discuss its future potential. 
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Chapter 2 

SR: Problem Formulation and 
Modelling 


In this chapter, we present a new approach toward solving the superresolution problem. 
Simplicity and a direct connection to the problem of single image restoration (from one 
observed image) are the main strengths of this approach. The various known methods to 
restore one image from one observation are easily generalized to the new problem of single 
image restoration from several measured images. We start our presentation with a new 
model of the problem and then turn to apply known restoration methods to the suggested 
model. The key to a comprehensive analysis of the classical superresolution problem is 
to formulate the problem and to model it as simply and as efficiently as possible. We 
start by presenting the problem to be solved and then turn to introduce an analytical model 
describing it. 


2.1 The Proposed Algorithm 

Our algorithm is applicable to spatial domain reconstruction, in which we have considered 
a space varying PSF. This algorithm can be used for reconstruction from different frames 
obtained by multiple cameras with different sensor element sizes. The recursive recon- 
struction of high resolution can be started from the very first frame and adding each new 
frame improves the resolution of the reconstructed image. We have used different type of 
blurring kernels for different frames in the simulation. 

After each recursion, we use our a-priori knowledge of the fact that each pixel value 
in the original image (the object) is nonnegative. Since a regularized least square solution 
doesn’t guarantee nonnegative pixel values, we simply force all negative pixel values to 
zero after each recursion, thus decreasing the MSE, This given an estimate that is closer 
to the object and acts as a better input for the next stage of recursion. Assume that after 
a recursion we have an estimate X {m, n) of original image X (m, n) and after making 
negative pixels values zero we get X' (m, n). Let MSSi and MSE^ be mean square errors 
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before and after using the a-priori knowledge respectively. We can write: 

MSEi = -^ E E [X (m, n) - X (m, n)] ' (2.1) 

m=ln=l 

where X {m, n) can take any positive or negative values. 

Using a-priori knowledge we can write; 

X' {m, n) = X (m, n) for X {m, n) >0 (2.2) 

= 0 f(yrX{m,n)<0 (2.3) 

I L L 2 

MSE2 = 72 E E [^' ^ «)] ^ ( 2 - 4 ) 

^ m=l 71=1 ^ 

2.2 Problem Formulation 

We are given N low resolution, blurred and noisy frames of different sizes Mk x Mk for 
1 < A: < AT. Problem is to reconstruct a high resolution image from these frames i.e. we 
have to find an optimal estimate of the original high resolution image. 

We now present a new approach towards superresolution reconstruction. We start with 
the problem model and go on to find a recursive estimate of the solution. Our main empha- 
sis shall be on how best to use the extra information obtained by adding each new frame 
in the recursion. There are at least two approaches to estimate the solution from the given 
number of constrained equations. The first approach [18] divides the given constrained 
equations into groups and makes an estimate by minimizing the sum of the mean square 
errors for each individual group. 

Since the larger the number of constrained equations, the lesser is the bias on the esti- 
mate, by the above approach, the estimation error or bias for each group will be more and 
any estimate obtained by minimizing their sum will give a more biased result. Our approach 
proposes to use all the constrained equations simultaneously to estimate the solution. In 
our algorithm, the dimensions of H and y always increase according to the dimensions of 
each new frame added as shown in (4.17) and (4.18). This is the basic difference between 
our algorithm and that proposed in [18]. 

In one method, we have applied the above algorithm without any type of denoising 
filtering. Since a small value of A (the regularization parameter, discussed later) gives 
good deblurring but poor denoising whereas a large A provides poor deblurring but good 
denoising ([13] and [6]), we have taken another approach, discussed in later chapters, the 
two variants of which both keep A low, and reduce the noise subsequently by subjecting the 
estimate to adaptive Wiener filtering [3] . The filtered estimate is used for the next recursion. 

2.2.1 A Discrete Formulation of Reconstruction 

We assume that an object at infinity is being imaged by N separate cameras after light from 
it has passed through a linear space varying medium. The medium is assumed to be locally 
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N 
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Frames 


Figure 2.1; Degradation model 

Figure shows the degradation model for the super resolution restoration problem 


space invariant, so that each input data frame captured is blurred only by a linear space 
invariant PSF, even though the PSFs that blur successive frames may be different. Further, 
additive noise is also assumed to have been introduced. 

The problem is to reconstruct a high-resolution image i.e. an estimate of the source 
image. Fig. 2.1 shows the overall situation of the problem under investigation, x is the 
object, from which we obtain N blurred ‘observation’ frames. These degraded frames are 
assumed to be obtained from the source image using the following equations: By the Fig. 
2 . 1 , 

00 00 

tk {i,3) =13 S {hf m,n) a: (m,n) ; 1 < A: < iV (2.5) 

m=— 00 n=—oo 

These blurred frames are shifted by (5**, 5^^) and the downsampled by {Zi^,Zjf). Finally 
we get N degraded frames y*. (i, j) by the following equations. 

Vk ihJ) ~ I'k {l^ik (2-6) 

oo oo 

yk(i,j)= h(i^zk+^zkJ^jk + ^jk’^>'^)*^(^>‘^) + ^k(iJ);^<k<N 

m="O0 n=— 00 

(2.7) 

where 6* is the blurring function for the kth frame, Z^^, Zj^ are downsampling factors, 
and 5,^, 5j^. are shifts in the i and j directions respectively, and x{i,j) is the source image. 
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2.2.2 Vector Matrix Formulation 

For Vector-Matrix notation we can rewrite, given N observation images {Yk}^^ , where 
each image is (in the general case) of a different size [Mj. x M^]. We assume that these 
images are different observations of the same high-resolution image X of size [L x L], 
where typically L > Mk fox 1 < k < N. More specifically, each observed image is the 
result of a shift, a linear space- variant blurring, and uniform rational decimation performed 
on the ideal high-resolution image X. We further assume that each of the observations is 
contaminated by nonhomogeneous additive Gaussian noise, uncorrelated across different 
observations. In order to treat the most general case, it is assumed that each observation 
is the result of different blur, noise, and decimation parameters. Translating the above 
description into an analytical model, we get: 

Using matrix notation, we rewrite (2.7) as follows: 


Yk = DkSkBkX -I- Bki 1 < k <N (2.8) 

where Bk is a [L^ x L^] matrix representing the linear space variant blurring operator ma- 
trix, Sk is a [L‘^ X L^] shift operator matrix, and Dk is an x matrix representing the 

decimation operator. In this model yk and x are lexicographically ordered. Ck is the addi- 
tive zero mean gaussian noise in the A:th observation with a positive definite autocorrelation 
matrix Wk“^ of size [M| x M|]. All these matrices (Bk,Sk,Dk,Wk,yk) are assumed to 
be known. Rewriting (2.8) we get. By grouping the N equations into one from the model 
equations described by (2.8) we get. 


yk = HkX + e (2.9) 

where Hk= DkSkBk. By grouping the AT equations into one from the model equations 
described by (2.9) we get, 

y = Hx-|-e (2.10) 

where 

y = [y? y? • • • yS]’^: 

H = [H? Hi' 
e = [e? ei ‘ 


In solving for x from (2.10) which is an inconsistent linear system of equations Ax = b 

where A is known to be an ill conditioned The least squares solution, x = (a'^ A^ A'^b, 
is not accurate enough due to ill conditioning of matrix A. In order to obtain a stable 
solution, a regularization approach is used which is discussed in the next chapter. 
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Chapter 3 

Inverse Problem and Regularization 


The desire to extract information regarding the structure of a signal or image given a noise 
corrupted, “blurred” version of the original is a common goal in many fields of engineering 
and the applied sciences. These inverse problems arise in fields as diverse as geophysical 
exploration, medical imaging, non-destructive testing, and radar signal processing. For ex- 
ample, a common signal and image processing problem is that of deconvolution where one 
observes a filtered version of a signal in additive noise and seeks to recover the uncorrapted 
original. 

While common enough in practice, problems such as these are notoriously difficult to 
solve. Most inverse problems are characterized by an unusually high sensitivity to pertur- 
bations in the data so that a small change in the measurements results in wild changes in 
the recovered signal. 

What is an Inverse Problem? 

Inverse Theory is concerned with the problem of making inferences about input data to 
a known system from observed output data (usually remotely sensed). Since nearly all data 
is subject to some uncertainty, these inferences are usually statistical. Further, since one 
can only record finitely many (noisy) data points and since physical systems are usually 
modelled by continuum equations that inherently imply bounded changes in output for 
bounded changes in input (at least geophysical ones are) no geophysical inverse problems 
are really uniquely solvable: if there is a single model that fits the data there will be an 
infinity of them. (A model is a parameterization of the system, usually a function.) Our 
goal then is to characterize the set of models that fit the data and satisfy our prejudices as 
well as other information. 

To make these inferences quantitative one must answer three fundamental questions. 
How accurately are the data known? How accurately can we model the response of the 
system? In other words, have we included all the physics in the model that contributes 
significantly to the data? Finally, what is known about the problem independent of the 
data? This is called a priori information and is essential since for any sufficiently fine 
parameterization of a system there will be unreasonable models that fit the data too. Prior 
information is the means by which we reject or down-weight unreasonable models. 

uIj jo 



3.1 Well conditioned systems 

The systems in which the mathematical problems fulfill Hadamard’s definition of well- 
posedness for which all the following properties hold: 


For all admissible data, a solution exists. (3.1) 

For all admissible data, the solution is unique. (3.2) 

The solution depends continuously on the data. (3.3) 

are well conditioned systems. 

The problems for which one or more of the above properties do not hold are ill posed 
problems. 

If a problem is ill-posed, one is usually not too much concerned with the violation 
of (3.1), although of course also existence of a solution (for exact data) is an important 
requirement. (3.1) can usually be enforced by relaxing the notion of a solution at least for 
exact data, while for perturbed data, the problem has to be “regularized” and hence changed 
anyway. 

Violation of (3.2) is considered to be much more serious. If a problem has several 
solutions, one either has to decide which one is of interest (e.g., the one with smallest 
norm, which is appropriate for some, but not all, applications) or one has to check the 
model for completeness and, if possible, feed in additional infonnation. It might happen 
that in a practical problem, the available data (even if measuured exactly at infinitely many 
points) simply do not determine the quantity which is sought. Then one has to “invent” 
additional measurements. 

The question of uniqueness is relevant in inverse problems where one looks for a cause 
for an observed effect; One is usually content or even happy with having a variety of pos- 
sible solutions (causes), since then, one can try to pick one which fulfills some additional 
criteria. 

Violation of (3.3) creates serious numerical problems: if one wants to approximate a 
problem whose solution does not depend continuously on the data by a “traditional” nu- 
merical method as one would use for a well-posed problem, then the numerical method 
becomes unstable. A (partial) remedy for this is the use of “regularization methods”, al- 
though one has to keep in mind that no mathematical trick can make an inherently unstable 
problem stable. All that a regularization method can do is to recover partial information 
about the solution - provide an inaccurate solution - as stably as possible. The “art” of ap- 
plying regularization methods lies in finding the right compromise between accuracy and 
stability. 
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3.2 The ill conditioned nature of restoration 

In the high resolution reconstruction problem where blur distortions are included, it is 
desirable to include the deblurring computation into the reconstruction process since the 
deblurring of input frames separately would introduce phase and high frequency distor- 
tions in the input frames, which is not desirable for high resolution reconstructions [13]. 
When the deblurring is incorporated into the high resolution restoration process, the overall 
reconstruction becomes unstable due to ill-conditioning of the system of equations. 

The problem of image restoration is the determination of the original object distribution 
X given the recorded images yj, and the point spread function 6. A common method of 
analysis of equations such as (2.7) is by operator theory. Given the space of functions that 
contains x and the space of functions containing the yk, the transformation or operator Tk 
that maps x into yk can be posed as 

Tkix}^yt,l<k<N (3.4) 

where obviously for images 

OO 00 

S S + + (3.5) 

m=-oo n=~oo 

By combining all the equations for 1 < fc < JV into one from (3.4) and (3.5); 

T {x} y, (3.6) 

The problem of image restoration is then to find the inverse transformation T~^ such that 

T-" {y} X (3.7) 

In a mathematical sense, the problem of image restoration corresponds to the existence and 
uniqueness of the inverse transformation. Both existence and uniqueness are important. If 
the inverse transformation does not exist, then there is no mathematical basis for asserting 
that X can be exactly recovered from y. (Although there may be no mathematical basis for 
exactly recovering x, there may be a practical basis for asserting that something very close 
to X can be recovered.) Problems for which there is no inverse, transformation, i.e., T“^ 
does not exist, are said to be singular. On the other hand , T"^(x) may exist but may not 
be unique; i.e. may not be a function. Finally, even if T~^ exists and is unique, it may 
be ill-conditioned, by which we mean that a trivial perturbation in y can produce nontrivial 
perturbations in x. That is, there exists e, which can be made arbitrarily smaU such that 

T"My + e} = x-f e' (3.8) 

where e' » e is not arbitrarily small and is not negligible. Thus an ill-conditioned prob- 
lem is one in which inherent data perturbations can result in undesirable effects in the 
solution by inverse transformation. As will be seen, image restoration problems are in this 
class at best and are frequently singular in addition. 
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3.2.1 Condition Number 


Consider the elements of the matrix H. If it is assumed that H is a non-trivial blur/shift/downsample 
operator, then it is reasonable to expect that the PSF represented by H is an image with con- 
siderable correlation between adjacent pixel rows and/or columns. Thus, the rows of the 
matrix H can be approximated by: 

h{i + l,j]m,n) = h{i,j\ m,n) + a [h {i,j‘,m,n) — h{i — l,j-,m,n)] . (3.9) 

which is to say, that in the matrix H, index iis a row index and (3.9) shows that the (i + l)th 
row of H is approximately a linear combination of the zth and (i — 1) rows. 

The relation (3.9) is of importance for the following reason: since the rows of the 
matrix H are approximately linear combinations of one another, the matrix H is very nearly 
singular. From matrix theory it is known that equations with near-singular matrices are 
extremely difficult to solve; i.e., they are ill-conditioned. This property of a matrix is is 
formally addressed in matrix theory by the condition number. Suppose that in the equation 

y = Hx, (3.10) 

if the matrix H is a square matrix, it is known from elementary matrix theory that the rank 
of H is equal to the number of nonzero eigenvalues of H. Let pmin and p^ax be the smallest 
and largest eigenvalues of H; condition number may be defined as 

c(H)= (3.11) 

\ Pm%n j 

From (3.11) it is obvious that as p^m. — > 0, the condition number approaches infinity. 

Thus, even though the matrix H may not be singular, near singularity is associated with 
very small eigenvalues and the solution of the resulting linear equations becomes difficult. 

In the presence of noise (whether sensor noise or computational noise) the matrix H can 
be singular within the bounds of uncertainty imposed by the noise. Solution of the digital 
restoration problem is thus tied to the solution of ill-conditioned (near-singular) systems of 
linear equations and computational methods become important. 

3.2.2 In Terms of SNR 

Concentrating on the former analysis, assume initially that the matrix H is square and 
nonsingular. An estimate of x can be generated by 

X = H'V, (3.12) 


and from (2.10) 


X = H"^ (Hx -1- e) 
= X -b H-^e. 
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(3.13) 

(3.14) 



Thus, the estimate is composed of two parts: the actual object distribution and a term 
involving the inverse acting on the noise. 

In an algebraic analysis, (2.10) is considered as a linear system with uncertainty or 
errors in the data; i.e., the data is “good” to a limited number of digits because of the 
measurement errors. If H has a large condition number, i.e., near singular, then the inverse 
will have very large entries since the determinant of H is small, and consequently the 
term H“^e will dominate the term containing the solution x, so that in the presence of an 
inherent error such as e there is a large inherent uncertainty in the solution. 

Now, assume that in (2.10) 


I y 11= II Hx + e ||< II Hx |1 H- 1| e || (3.15) 


where |1 y 1| is the Euclidean norm of vector y. Now assume that H preserves energy in x; 

II Hx 11=11 X II, (3.16) 

SO that the relative signal-to-noise ratio (SNR) (defined in terms of norms) is 


Now in the solution. 



(3.17) 


I X 11=11 x + H II < II X II + II H ^ nil e |1, (3.18) 


where an Euclidean norm for both vectors and matrices is used. If the norms of the two 
terms in (3.18) are compared. 


I ^ I 

H-i ||||e|' 




(3.19) 


even though H preserves energy in x, does not demonstrate such properties in general. 
Indeed, it is typical of an ill-conditioned problem that 


II 11» 1. 


(3.20) 


Thus, even though the original data’s signal-to-noise ratio was a, the solution possesses a 
signal-to-noise ratio of ^ and ratio of these two gives the deterioration in signal-to-noise 
between data and solution: 

Deterioration factor in SNR 
from data to solution 

For example, suppose that originally || x ||= 100, || e ||= 1, and || H“^ ||= 10. Then 
os= 100, ^0= 10, and deterioration factor = 10, Thus, even very good signal-to-noise ratios 
can be deteriorated by ill-conditioning. (|| ||= 10 would be considered very mild in a 

restoration problem; || ||= 1000 or more is typical, meaning the noise would dominate 

the solution for the hypothetical case posed in this example.) 


) = ^ =11 H-' II . (3.21) 
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3.3 Regularization 

As we have seen, the equation Tx = y is generally ill-posed when T is a linear operator. In 
particular, (2.7) is ill-posed when the kernel h {., .) is square-integrable. If we are to obtain 
useful approximate solutions to either of these equations, we must modify the problem in 
such a way that the solution becomes less sensitive to small perturbations in the data. At 
the same time, the solution to the modified problem must be close to the solution of the 
original problem. The process of modifying the original problem to achieve an acceptable 
compromise between these two conflicting goals is referred to as a regularization. We 
briefly discuss the principles of regularization in this section. 

The ill-posed nature of the equation Tx = y is a consequence not only of the nature of 
the operator T but also of the domain I?(T) over which T is defined. Hence, we may be 
able to regularize the problem at hand either by modifying D(T) or by modifying T itself. 
Regularization involves replacing the operator T by an operator T. Suppose further that 
“true” y is perturbed, yielding an observed signal y satisfying 

I! y - y ll< £. (3.22) 

We seek a least-squares solution to the modified equation 

tx= y. (3.23) 

If the modified operator T has a bounded pseudoinverse T \ the least-squares solution 

X = T+y (3.24) 

will be stable and is a potential candidate as an approximate solution to the original equa- 
tion. We seek modifications of the original operator T that yield approximate solutions of 
the form (3.24) having an appropriate balance of accuracy and stability. 

We express the original least-squares problem in the form 

Tx ~ y. (3.25) 

as it stands, this least-squares problem still suffers firom instability and must be regularized 
if a useful approximate solution is to be found, A regularizer for this problem is defined as 
follows. 

Definition: A regularizer for (3.25) is a one-parameter family of linear operators 
{Ry, z/ € A} , where A is a set of positive real numbers with 0 G A, satisfying 

(а) For each z/ > 0, is bounded, and 

(б) lim^_+o 11 Ri/Tx - x 11= 0 for x in the range space of 

Thus, R,, is a bounded approximation to T^ and can therefore be used to obtain a stable 
approximate solution to (3,25) with an appropriate choice of the regularizing parameter v. 
For a given value of u, let us take 

Xi, = Ri,y (3.26) 
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as the approximate solution to (3.25). We 


then have 


li X, ~ X |j=|( 11^5^ _ Tty -f R,y-R^y 11 

^11 H^y - R»,y II + (I R^y - Tty || 
<11 R^ lle+ll R^y-Tty I) . 


The second term is thp ^ 

approximation ,o Tt). and itT" t' «> the fact that R„ is only an 

of fidelity due to the applicatiorof '>'“ff«ave loss 

of the observation error. This term dn ^ ^ ii ^ ^ ^ amplified version 

always blows up a fact that -n u not generally vanish; indeed, as i/ — > o it almost 

stability (controlled noise atnDlmS-™rT^*T^” rogulaiization error) and 

poor fidelity, while too small a v^e"!! ' d° ““»■> of 

ill effects of noise amplificatiol'"^''^ insufficient regularization to overcome the 

tem!f h“Sr eqntta'^thttS“''' “>0 ill-conditioned sys- 

Iteration, the number of eouations restfianzation function. In the form of 

rated for improved Jutil n h “ ? ’ ” “fo™ation can be incoipo- 

an improved solution bv uHlv otems^es shown in Hg.i, it is desired to^et 
conditioning problems For thT® ^“““^^.““snremeuts while overcoming the Ul- 

rogularization rosulJiraZ^ 'P’”^’ ‘f’' ™ effective 

aro incoTiorated “fonnatron is recurstvely updated while new measurements 


y Squares Regularization for Recursive Reconstruc- 

?on to mi”!!!?!.! legulaiization approach 


is used which finds the solu- 


R.(x) = ||Hx:-yf -(- tW (3 27 ) 

taortox^^e! (2.10). If we 

7(x) = A||x-ci|%h::AT""‘^^‘°’" 

ow the overall objective functional to be minimized becomes 

R(x) = ||Hx - y f + A||x - cjp V 3 _ 2 g) 

and the solution x that minimizes R(x} is 

X = (h^'H + AI) "" (H'^y + Ac) . ^ 29) 
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where A is a regularization parameter (same as u in previous section) which is adjusted to 
get a satisfactory compromise between stability and accuracy of the solution. A controls 
the effectiveness of deblurring, high resolution reconstruction as well as noise reduction as 
discussed in previous section. 
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Chapter 4 

SR: Recursive Reconstruction Procedure 


In the first section of this chapter, we discuss adaptive Wiener filtering and in the second 
section, we present the recursive reconstruction procedure to be applied for getting a high 
resolution image. From ( 3 . 29 ), to get an estimate of x we must invert ^H'^H + AI^ . Here 
for implementation, we have considered a small dimension image, so that the inversion 
in ( 3 . 29 ) does not take too long to compute. The computation would be quite heavy for 
larger dimension images. In such cases, we could use iterative methods such as the steepest 
descent (SD), conjugate gradient (CG) [ 15 ] etc. Bellman [ 4 ] describes an iterative approach 
to the solution of the ill-conditioned system of linear equations with an iterative update of 
the regularization function. 

4.1 Adaptive Wiener Filtering 

Wiener filter is a minimum mean square estimator. Let /(ni, 712) and q{ni, 112) be arbitrary 
random sequences, where q{ni,n^ is a noisy version of /(ni, 712) related as follows: 

5(711,712) = f{ni,n 2 ) + v{ni,n 2 ) ( 4 . 1 ) 

it is desired to obtain an estimate, /( 7 ii, 772 ) of / (711,772) from 5(7x1,722) such that mean 
square error = £^{[/( 72 i, 7x2) — /(txi, 722)]^} is minimized. 

We use adaptive Wiener filtering because simple Wiener filtering blurs the image exces- 
sively. When an adaptive image processing method is applied to the problem of restoring an 
image degraded by additive random noise, it is possible to reduce background noise with- 
out significant image blurring. Two approaches to adaptive image processing have been 
developed. In one approach, called pixel-by-pixel processing, the processing is adapted at 
each pixel. At each pixel, the processing method is determined based on the local charac- 
teristics of the image, degradation, and any other relevant information in the neighbourhood 
region centered around the pixel. Since each pixel is processed differently, this approach is 
highly adaptive and does not suffer from artificial intensity discontinuities in the processed 
image. In another approach, called subimage-by-subimage or block-by-block processing. 


19 




nif+m^ 


ffif 


Figure 4 . 1 : Wiener filtering 

Noncausal Wiener filter for linear minimum mean square error estimation of f{n\, 712) from 
g(7ii,n2) = /(ni,n2) + u(ni,Ti2). 


an image is divided into many subimages, and each subimage is processed separately and 
then combined with the others. 

Most adaptive restoration algorithms for reducing additive noise in an image can be 
represented by the system in Fig. 4 . 2 . From the degraded image and prior knowledge, 
some measure of the local details of noise-ffee image is determined. One such measure is 
the local variance. A space variant filter g{ni, 712) which is a function of the local image 
details and of additional prior knowledge is then determined. 

The space-variant filter is then applied to the degraded image at the local region from 
which the space- variant filter was designed. When the noise is wide-band, the space- variant 
ginifti-z) is low pass in character. In low-detail image regions such as uniform intensity 
regions, where noise is more visible than in high details regions; a large amount (low cutoff 
frequency) of lowpass filtering is performed to reduce as much noise as possible. Since little 
signal variation is present in low detail regions, even a large amount of low pass filtering 
does not significantly affect the signal component. In high-detail image regions such as 
edges, where a large signal component is present, only a small amount of low pass filtering 
is performed so as not to distort (blur) the signal component. This does not reduce much 
noise, but the same noise is less visible in the high-detail than in the low-detail regions. 

A number of different algorithms can be developed, depending on which specific mea- 
sure is used to represent local image details, how the space variant g{ni,n2) is determined 
as a function of the local image details, and what prior knowledge is available. One of 
the many possibilities is to adaptively design and implement the Wiener filter. As Fig. 4.1 
shows, the Wiener filter requires knowledge of the signal mean m/, noise mean my, signal 
power spectrum Pf{u}i,u}2), and noise power spectrum Instead of assuming 

a fixed m/, rriy, Pf{oJi,uj2) and P/(wi,W2) for the entire image, they can be estimated 
locally. This approach will result in a space-variant Wiener filter. 

We first assume that the additive noise u( 7 ii, nz) is zero mean and white with variance 
of al- Its power spectrum is Py (wi W2) then given by 

(Wl,W2) = (4-^) 

Consider a small local region in which the signal f{ni,n2) is assumed stationary. Within 
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Figure 4.2: Adaptive filtering 

Topical adaptive image restoration system for additive noise reduction.) 


the local region, the signal /(ni , 712) is modelled by 


/ (ni , n2) =mf + (TfU) (ni, 712) 


( 4 . 3 ) 


where rn / and cr/ are the local mean and standard deviation of / (tii , 712) , and w (tii , 7^2) is 
zero-mean white noise with unit variance. 

In (4.3), the signal / {nunf) is modeled as a sum of a space- variant local mean ruf 
and white noise with a space-variant local variance a'j. Within the local region, then, the 
Wiener filter G {u}i,wf) is given by 


^ i^/(wi,W2) _ 

“ F/ (a;i, 072) + F« (071,072) aj + at 


(4.4) 


from (4.4), g{nun2) is a scale impulse response given by 


p (711,712) 



5(711,712). 


( 4 . 5 ) 


From (4.5) and Fig. 4.1, the processed image ^(711,712) within the local region can be 
expressed as 


p(ni,7Z2) = TTT/ + (9(7^1,712) 7 7n/) * 


771/ -t- 




0^f+(^v 


(9(771,772) - 777 /) 


( 4 . 6 ) 

( 4 . 7 ) 
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If wc assume that m/ and aj are updated at each pixel, 

2 / 's 

Pin,. .,) = m,(n„n,) + n,) '''■«) 

Noting that ruf is identical to niq when is zero (and it is indeed zero by assumption,), 
we can estimate m/(ni , n^) in (4.8) from q{ni, 7^2) by 


mf (ni,n2) = 


m+M n^+M 

(2M + 1? ^ ^ 

-t 1) k2=n2-M 


(4.9) 


where (2M + l)'"^ is the number of pixels in the local region used in the estimation. 
Since = a j + oj may be estimated from g{ni,n 2 ) by 




(ni, na) - <^, if (ni, ^2) > crl 


0 


otherwise 


where 


^qinuTiz) = 


ni’^M n2+M 


S Y, (9 {ku M - m/ (ni, na)) 


(2M + 1) ki=rn-M ki=n2-M 


(4.10) 


(4.11) 


4.2 Recursive Reconstruction Procedure 

In this section wc show how to make recursive the procedure to update the estimate. Since 
the solution given by (3.29) docs not guarantee a reconstruction with nonnegative pixel 
values, wc use the a-priori information that x(z) > 0 for I < i < L^, and, after each 
recursion, set all negative pixel values to zero thereby decreasing the MSB. This obviously 
gives an estimate that is closer to the desired one and acts as a better input for the next stage 
of recursion. The Wiener filter output for the estimate Xk is denoted by Wk. 
Initialization : A=0.01. 


c(l) = 0; 

(4.12) 


(4.13) 

y (1) = yi; 

(4.14) 

X, = (tf (1)’' H' (1) + Al)"' (h' (!)■'/ (1) + Ac (1)) . 

(4.15) 

Update: 1 < A: < 16 


C (k) = Wk-l; 

(4.16) 

H'(k) = [H"'(k-l) Hj]''; 

(4.17) 

/ W = [y™ (k - 1) yf]''; 

(4.18) 

Xk = (h' (k)’' H' (k) + Al)"' (k) y (k) + Ac (k)) . 

(4.19) 
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Figure 4.3: Recursive reconstruction 

Recursive reconstruction model for a high resolution image from N degraded undersampled 
frames. 


In the next chapter, we discuss the implementation and compare the results obtained by 
different algorithms. 
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Chapter 5 

Implementation and Discussion 


In this chapter wc discuss about two algorithm. First one is for reconstruction of a high 
resolution image from several degraded frames, which is much computationally expensive. 
Second algorithm is for rcconstiuction of a high resolution image from a degraded frame 
blurred by a separable point spread function. 


5.1 First Algorithm: Computer Simulation and Results 

For simulation on MATLAB wc use ‘cameraman.tif’ image of size 32 x 32. The reason for 
taking a small image is to Ics.sen the simulation run time. From this picture we obtain 16 
undcrsamplcd, blurred and noisy frames of size 8x8. The blurring kernel Hk; 1 < fc < 16 
for each frame is different. 

Hence the computer simulations use 1 6 frames of 8 x 8 noisy, blurred and undersampled 
images as the input frames in order to reconstruct 32 x 32 high resolution image. The low 
resolution blurred input frames arc obtained as follows. First 32 x 32 original image is 
blurred by different types of 2-D Atmospheric turbulence blur and Motion blur. Then the 
image is downsampled with various sampling references to get low-resolution frames with 
different shifts. Then Gaussian distributed noise are added in each frame. Fig. shows the 
16 input frames of size 8x8. For simulation we have taken the value of regularization 

parameter A as 0-01 . i, uu 

Wc show the improvement in the estimated image after adding each new frame both by 

visual comparison and by plotting the mean square error. , , . r 

We compare the results from our algorithm with those obtained by [18] on the basis of 
mean square error of the estimate. We have given equal weight to each past frame: the use 
of a weighted least squares approach and/or a smoothing operator [18] might improve the 

3r0Si^ltiS 

For purposes of comparison, we consider various different cases as follows: 

Case 1. SR images obtained without Wiener filtering of the estimate obtained from 

(4.19): then (4.16) becomes c (k) = Xk_i. 
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Case 2. SR images obtained using (4.16) to (4.19). But now, Wiener filtering is applied 
only to prepare the input data for the next incursion, while the partially re-solved image 
after the kth recursion is taken to be Xfc, 

Case 3. SR images obtained using (4.16) to (4.19). This means Wiener filtering is 
applied to prepare the input data for the next recursion, and the partially resolved image 
after the kth recursion is taken to be . 

Case 4. SR image obtained by method suggested in [13]. 

Case 5. SR images obtained by using the recursive algorithm given in [ 18]. 



(c\ 

Figure 5.1 : (a) Original image of sizes 32 x 32, (b) Input degraded frames of sizes 8x8, 
(c) Reconstructed image according to case4. 

In results we have discussed two observations. In first one the reconstruction of high 
resolution image is started from the first frame. In second observation we have started 

reconstruction assuming six frames available. . . , . 

Fig, 5.1 is common for both the two observations. Fig. 5.1(a) shows the ongmal image 
of size 32 x 32. Fig. 5.1(b) shows the input undersampled, blurred and noisy frames of size 
8x8. Fig. 5.1(c) shows the constructed image described as in case4. 


5.1,1 First Observation: Computer Simulation and Results 

The result for case4 is same irrespective of input MSEs order. Now we^scuss abouUhe 
observation when the MSEs of the input frames are m random fashion. Case2 is still best 
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on the buses ot coinpuring the MSE among all the cases. But the output in case2 is some 
blurred due to single time denoising by the Wiener filtering in each recursion. 



i'igure 5.2: Random ordered Input MSE: MSE comparisons graph for different cases. 


5.1,2 Second Observation: Computer Simulation and Results 

Here wc discuss about second observation. It is observed that if the recursion is com- 
menced, not from frame k = 1, but at a frame ko further down, and the proposed algorithm 
(casc2) is modified as described below, the MSE at k = 16 is lower than that seen when the 
same algorithm is applied from k = 1. 

Initiali'/alion: A={).0L 


c (ko) = 0; 

(5.1) 

H'(ko) = [Hr 

(5.2) 

y'(ko) = [y? y? ••• yS]"^; 

(5.3) 

Xko = (H' (ko)*^ H' (ko) + (H' (ko)^ y' (ko) + Ac (ko)) . 

(5.4) 

Update: k© < k < 16 


c(k) = Wk_i; 

(5.5) 



(5.6) 

(5.7) 

(5.8) 


H'(k)=[H'(k-l)^ 
y' (k) = [y (k - 1 )^ yT]^; 

X, = (H' (k)"^ H' (k) + Al) (h' (k)^ / (k) + Ac (k)) . 

'Fhis is demonstrated in Fig. 5.4 and Fig. 5.5 for a starting frame of k = 6. The reason 
for this phenomenon is iwt conclusively evident, but it is likely that, after adding a certain 
number of frames, noise amplification starts in the partially resolved SR image [13] and 
this noise propagates through the recursion. The curve for case2 in Fig. 5.2 after k = 10 
suppt>rts this hypothesis. 

At pivsent. further study of this phenomenon is required before a clear picture emerges, 
and (.{uesiitms such as whether there exists an ‘optimal’ start frame number, can be an- 
swered. 


5.2 Second Algorithm: Computer Simulation and Results 

Now further we discuss about a new computationally efficient method for separable PSF 
taking the fad that 2-D shift matrix and 2-D downsampling matrix are always separable 
and for getting a stable and unique solution of the ill-conditioned equations we use empir- 
ical rcgulari/ed least square mcthtxJ. The direct invenion of + AI^ raises several 
difiicuities which render the solution of (3.29) impractical. The main one is the computa- 
tional burden required here. In a typical video application the size of the matrix is x 
(c.g. for L. -HMX). wc will be required to invert a 10® x 10® matrix). Since computational 
complexity for finding the inverse of a matrix of size [N x N] is of Af® order. Eq. (3.29) 
is having matrix of si/c x which will have computational complexity of I® order. 
For overcoming such situation, in the proposed algorithm we have exploited the separable 
property of shift matrix S and downsample matrix D. hence if the blurring matrix B is also 
separable according to following equations 
D = Di 0 Da 
S = Si 0 Sa 
B = Bi 0 B 3 

Using above relations rewriting (2.8), 

y = DiSiBiXBaSaDa + E (5.9) 

where each entity in above equation is a matrix. Bi, Ba, Si and Sa each is of size [L x L] 
and D, , D J arc of .size (M x L] and Y, E are of size [M x M]. Rewriting (5.9) we get, 

Y = HiXHa + E (5.10) 

where Hi = DiSiBi and Ha =* DaSaBa- If we use regularized least square method for 
getting estimate of X from (5.10), wc getthe sameestimate equation as (3.29), due to which 
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m,r m..livc of a.nsWcring PSF to be scpiirable seems useless. Hence for solution 

possible and compuumonally cllicicnt, wc propose an empirical formula for separable case 

which is inspired hy (3.29). 

I-rom H - Hi Ha for gelling empirical solution, a partial estimate Xi of X can be 
given as follows considering only Hi, 

X. = {H?H, + Al)-*(HrY). (5.11) 

similarly t\>r considering only Ha another partial estimate Xa of X can be given as follows 

Xa = (YHJ) (H^Ha + Al)"^ , (5.12) 

Now considering both Hj and Ha simuUancously and combining above both partial esti- 
mates. wc can write a empirical estimate of source image as follows: 

X = (HjH, + Al)"' (HfYI^) (h^H, + Al)'* (5.13) 

Now for obtaining estimate from (5. 13), we have to invert two matrices each of size LxL. 
That is easily practically implcmcntable. 

In next section we will compare the results obtained by actual regularized least square 
method and proposed empirical least square regularized method. 

In this section wc have compared the results of the proposed algorithm with two differ- 
ent cases. In first case the estimate of the high resolution image is found by interpolation, 
deshifting. deblurring and further denoising by Wiener filtering or median filtering. In the 
second case the estimate of the original image is found by actual regularized least square 
method given in (3.29). which is quite computationally expensive. 

For simulation on MA'HLAB wc use ‘caraeraman.tif’. We have taken the size of image 
as 256 X 256 for comparison with first case and 32 x 32 for comparison with second case. For 
getting degraded frame first wc blur the original image by 2-D motion blurring, then shift 
it (here for simulation wc have taken zero shift in both directions.), then downsample it by 
2 and finally add white Gaussian noise. We compare the estimate of the first case with the 
estimate of the proposed algorithm for two different values of noise variance[cr^ = 0.001 
and «? = 0.0 1 1 . The estimate by the first case deteriorates very fast with the increasing noise 
variance which is also clear from Fig. 5.6c, Fig.5.6d, and Fig. 5.7b. While the estimate by 
proposed algorithm remains consistent as shown in Fig. 5.6e and Fig. 5.7c. 

For second comparison wc use image of size 32 x 32. The reason for taking small size 
is computations in second case. By Fig. 5.8 it is clear that the proposed algorithm is equally 
good as case second with drastic change in computational time. 

The simulation run time for 32 x 32 image is more than 30 minutes by (3.29) and 
restoration of same image by (5.13) takes few seconds. When the image of size 256 x 256 
is restored by (5.13), the simulation run time is even less than 1 minute, which is really 
remarkable. The restoration of the same image by (3.29) is unimaginable. 
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feua 5.5. Recursive reconstruction : MSE comparisons graph for different cases. 





Figure 5.6: (a) Original image of size 256 x 256, (b) Input frame of size 128 x 128 with 
noise variance a^^.OOl . (c) Wiener filtered reconstructed image for first case, (d) Me- 
dian fiUered reconstructed image for first case, (e) Reconstruction according to proposed 

algorithm . 
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Figure 5.7: (a) Inpul frame of sj/« 128 x 128 with noise variance cr^=0.01, (b) Median fil- 
tered reconstructed image for first case, (c) Reconstruction according to proposed algorithm 



Figure 5.8: Comparison in second case and proposed algorithm:(a) Onginal image of size 
32 X 32, (b) Input degraded frame of size 16 x 16, (c) Estimated image by eq. ( . ) ( ) 
Estimated image by cq. (5.13) 
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Chapter 6 
Conclusion 


in this chapter w.c conclude the thesis and discuss about the future work which can be 
priK'Ccded turther. 

'ITtc pniptjsed first algorithm gives better results but its computational complexity is 
substantial. ’Ihis catt perhaps he reduced by using some iterative methods such as SD and 
CG etc. 'I’he LMS algorithm for fast reconstruction given in [18] will give better results 
if we use our proposed recursive algorithm. In the second proposed algorithm we have 
restored an image trom a single degraded frame blurred by separable point spread function 
which is very much ctimputationaily efficient. 


6J Future work 

As discussed in chupler5, if we start recursion from other frame but not the first frame, the 
optimum start frame has to be find out by keeping compromise between deblurring and 
denoising in the final reconstructed SR image. Second work can be extended to get a high 
rcsolutitm image from several degraded, undersarapled frames for separable point spread 
function which would be a good achievement as computational complexity point of view 
for applicalit>ns in astronomy. 
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